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PROGRAM FOR THE DESIGN OF AN AXIAL COMPRESSOR STAGE 
BASED ON THE RADIAL EQUILIBRIUM EQUATIONS 



Kyriacos D. Papailiou 



1. Introduction 

The present report describes a program which has been established 
for the design of an axial flow compressor intermediate stage using the 
radial equilibrium equations as described in ref. 1 (pages 439 - 454). 

An intermediate stage is defined as one for which inlet and outlet 
velocities are essentially the same. 

The flow is assumed compressible and axisymmetric. The inlet 
conditions, some of which are data to the program, are assumed to have 
been produced in preceding stages and are referred to ambient condi- 
tions). Energy gradients at the inlet, which may be caused by 
non-uniform energy addition in previous stages, as well as entropy gradients 
at the inlet, which may be due to non-uniform energy dissipation in 
previous stages, are taken into account. 

The losses and the corresponding entropy increase through the rotor 
and the stator are calculated according to reference 3, and non-uniform 
energy increase in the rotor can be specified. 

The inner and outer walls may have arbitrary shape. However, when 
the curvature of the streamlines (in the way described in reference 1.) 
is taken into account, then the inner and outer walls, although maybe 
tapered, have to be straight. 

The problem is considered from the designer's point of view. Con- 
sequently, specified will be the work distribution and the quantity 
v ui +Vu2 

will be specified along the radius (see symbol table for symbols). 

The quantity ^ u l + ^ u 2 gives the theoretical reaction factor for an inter- 
u l +u 2 

mediate stage with constant axial velocity for incompressible flow. This 
quantity was chosen instead of the actual reaction factor as the expression 
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of the actual reaction factor is rather involved in the general case 
(see Appendix A). The thermodynamic diagram of the process through 
the stage is given in Figure (l) while the general layout is given 
in Figure (2) . 

In the following the theory will first be developed and then 
the program will be described. Symbols are defined in the symbol 
table and the FORTRAN Symbol Table . 
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Ef<r. 2. GE^LRni LAYOUT 

2. General Equations Used (See Ref. l) 




dV a 2 2 

dir- + FV a + G = ° 



( 1 ) 



where 



F = - 



_2 3Vr _1 1 dS 

V & dz Cp cos^A. dR 



G = 2 



S ( R V _ 2 SH , 1 






2x dS 



( 2 ) 

(3) 



with 



2 3V r _ 2K m _ + 6r ('ll 

V a dz cos3A 

The minus sign holds for stations (e) and (?) and the plus 
sign for station (?)* 



*Note that there is a discrepancy between the way equation (4) is inter- 
preted here and in ref. 1. However, positive curvature is considered here 
in accordance with the way equations (l), (2) and (3) have been developed. 
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(5) 



where tan \ e = tan = tan X 2 = _ 

Rq + L r 

L = — (See Fig. 2) ( 6 ) 

2 

K takes the values 4 to 6 . 

The solution of equation (l) is given as follows 



V 



2 



a 





R R 

JFdR JFdR 

B R h Rh 

e - e 



J g e 

R h 



R' 

J FdR ’ 

% 



dR 



(7) 



3. Dimensional and Non-Dimensional Quantities 



The basic quantities used are expressed as follows 
Lengths in ( ft) 
velocities in (ft/s) 
angular velocity u) in ( rad/ s ) 

enthalpies in (ft 2 /s 2 ) or 

. ,ft - lb N , ft 2 N 
entropies in (- 55 ^ 35 ) or (^or) 

the specific heat Cp in or 

the density in (slug/ft^) 

Consequently F has dimensions (ft“l) 

G has dimensions (ft/s 2 ) 

Calculations will be formed with non-dimensional quantities. As 
reference quantities will be used: The angular speed cd, the mean radius 

Rjjp where 

% + Rh-. 

Rm = -S (7a) 
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the atmospheric pressure, temperature and density. We shall call all 
non-dimensional quantities starred quantities and denote them with a 
star. 

- Lengths will be non-dimensional over P^. Consequently 





- Velocities will be non-dimensionalized over Consequently 

(where V ref = u;!^) 



V a 

V * = — - 
a a>R m 



W, 



W * = 
u 



u 



ujR, 



L m 



U* 



U 

m 



R* 



- Enthalpies will be non-dimensionalized over <ju c 




h* 



h 

A > 2 



H* = 



H 



2 2 

At 



- Entropies will be non-dimensionalized over C^ 

- Densities will be non-dimensionalized over the 

- Pressures will be non-dimensionalized over the 



atmospheric density 
atmospheric pressure 



4. Problem Formulation 

Using non-dimensional quantities (in the described way), we arrive 
at the following equations 



«( V a *> 2 

dR* 



+ F* (V *) + G* 



0 



(la) 



where 



_ 2 dV r * 
V a* Sz* 



1 ^(S/Cp) 

cos^l SR* 



(2a) 
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(3a) 



G* 



Vu* S(R*V„*) 
p* SR 



2dH* 

SR* 



+ (2H* - (V u *f) 



S(S/Cn) 

SR* 



and 



F* = FR m 



G* = G/o^R m 



<• 

a 



(8) 



We shall consider now Figure 3. 



Considering stations 




© 




Fig unL -3. 



we can see that the angle X can be found from 



tanX 



AR2 R3 - Rl 
2L 2L 



and still the curvature be calculated using &r at station (2) which is 
now negative measured from the straight line (T) - (T), instead from the 
straight line (J) - (^)- The curvature expression remains the same if 
this sign change is taken into account. 
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5 • Necessary Inlet Conditions 



Quantities will be specified at equal radial distances at the inlet 
The distribution of total enthalpy H and the distribution of entropy s 
must be specified at the inlet. This distribution of entropy and total 
enthalpy will be considered to have been developed through an adiabatic 
process (previously existing stages), starting from uniform atmospheric 
conditions. The atmospheric total pressure and temperature must be also 
given. Entropies will be measured considering as origin the atmospheric 



entropy . 




f~~/ (j-. 

Ut/TV PXL/C 7/orf. 



Considering an adiabatic compression 
from (a) (atmospheric conditions) to 
® (conditions at the stage inlet) 
along a stream surface, we have 
(see Fig. 4 for symbols) 

S B - S B* = C p ln <TT.> = V S A = &S 

-D 

or 

H b * = H b e‘ AS / C P 



from where Hg* can be calculated. For an isentropic process we have 
also (from state A to state B*) 

P tB* P tB H B * v/tY-i) / H b e- AS / C p \ T^Tj 

Pt A ' Pt A = ( « A = V S Tt A / 

from where the total pressure at the compressor inlet can be calculated. 
The stagnation density is given then as 

pt _ c p p t 
t Rg T t R g H t 
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6. General Method of Solution 



As can be seen an iterative process is necessary to solve the 
problem. A first approximation of the streamline position is assumed 
and on this basis the radial distribution of the axial velocity V * 

a. 

is calculated at stations (l), (2) and (3), that satisfies the 
continuity equation on the whole. 

Then, a new approximation of the streamline position is achieved 
by requiring that the same mass flow passes through individual stream- 
tubes at stations (2) and (3), the radial positions at station (l) 
remaining the same. 

A new V & *- distribution is then calculated and the iteration continues 
until the specified number of iterations is achieved. The printed 
error of each iteration gives us an indication of the convergence of the 
procedure. 

Having established the equations to be solved, a general flow diagram 
is given in Table I describing the general layout of the program. 

7. Detailed Calculations (Following flow diagram of Table I). 




Block Q 




equidistant radii are considered 



Ft<*. J. 



Block Cg ) 



Block 



The total enthalpy increase inside the rotor is 



AH = u> (Bg V U2 - Rp V U1 ) 



or AH* = 



AH 
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The reaction factor is 



V * V U2 V ul /«R m + V U2 /«R m V Ul * + Vu2 * 

EF = ~V r °2 V™m + U 2 / “ V + V 

or finally 

RF(R * + R *) - AH*/R p * 

v * = - £ _ 

U 1 1 + R */ R 2 * 

RF(R * + Rp*) + Z\H*/R * 

V * = - — _ 

U 2 1 + R */ R x * 



(9) 

( 10 ) 



V * = V 
V U3 v u]_ 



* 



Block G ) 

Calculations are performed vathout losses in the stage nor curvature 
effects for the first iteration. Corresponding terms are set here 
equal to zero. 



Block® 

The curvature term is calculated according to equation (4) rewritten 



here 

CRTERM 2 Sv r + &R 

Rm V a 9z l 2 



( 11 ) 



Note that in the way Vavra describes the curvature in page 453 of ref. (l), 
there exists an inconsistency in sign with the derived equations of motion. 
In fact the plus sign belongs to station z = 0 of ref. (l) or stations (l) 
and (3) here, while the minus sign to station (2) here (if is taken to 
be a positive quantity) . The expression for the curvature term then becomes 
2 dVj. §r 

CRTERM = - — ^4 R = + 2K R ra = 

V a oz m - l2 




In this expression 6 r* is positive in the way shown in fig. 3 and (+) sign 
applies to stations (l) and (3) while (-) applies to station (2). 
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The losses in total pressure are calculated in a subroutine, the 
theoretical basis of which is given in appendix B. This was decided in 
order to have the freedom of introducing any loss- correlation we desire. 

The subroutine will furnish to us the decrease in total pressure non- 
dimensionalized over the atmospheric pressure. 

The resulting entropy increase will then be calculated as 




where P Q corresponds to the total pressure level without losses. For the 
rotor this total pressure is P-tl' ( see fig* l) which is the corresponding 
total pressure that would result if the addition of enthalpy H2 - were 
done isentropically. Then we have 



jV 

p tl' 



% 

Hi 



V 

Y-l 



or 





Y 

T^i) 



and 



S 




1 





(13) 




For the stator P Q comes to be equivalent to P^ 2 and the increase in 
entropy from station (2) to station (3) is given as 

p t 2 - (APJ, 



AS -^g . 
jr- = - - 7 * In 

°p °p 



' o S 



t 2 



or 



S 3- S 2 



= s 3* " 



S 2* 



= _ M 






In 



S 3* 



S 2 - 



(15) 
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Block 



& 



The subroutine DERIR (see description in ref. (4)) is used to find 

dH* oS* 



the derivatives of a function given at discrete points. Thus 
d(R*V u *) 



and 



dR* 



are calculated. 



dR* ’ dR* 

Then the functions F* and G* are calcu- 



lated as given in equations (2a) and (3a). 

Block(j5) 

The integrals of equation (7) are calculated (FUNC1 and FUNC2) and 

y P 

the distribution of (v a ) is considered in the following form 



(V *) 2 = FUNC1 ( V„, *) 2 + FUNC2 

a a h 



An iteration is initiated using as starting V a ^* val ue one which does not 
make the quantity (V a *)^ negative. The iteration scheme is described in 
ref. 4 (iteration subroutine ITERN). 

For each value of (V^*), the (V a *) - distribution is calculated 
for the station under consideration. The iteration process ends when the 
continuity equation is satisfied. 

Block Q 

The continuity equation is considered in the following way. Consider 
the actual axial velocity distribution at a section, V a * = V a *(R*). Then 
the corresponding mass flow is 



Q* 



,-.R+ 



real 



Rh 



2-np*V a *R*dR* 



(16) 



Consider now a hypothetical situation where the part of the inviscid flow 
distribution of V * is extended inside the wall boundary layer regions 

a. 



up to the walls 
where 



Then using this distribution we get a mass flow rate Q*, 



Rt 

Q* = J 2rrp*V R*dR* 

R^ a inviscid 



then 



Q* > Q*. 



real 



and ^ai nv j_ sc i(j ( R* ) differ from V a *(R*) only inside the wall boundary layer 
regions. 
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We can say now that normally a flow rate Q* would pass through the 
area given,if it were not for the wall boundary layer presence. Defining 
a blockage factor K b as 



Q*real 

Q* 



( 17 ) 



then 



^*real 



A 



K1 



= Q* = | 2np*V * 

J Rv a 



R* dR* 



(18) 



inviscid 



where the mass flow rate has been non-dimensionalized as follows 



Q* = 



Q 






In the program the inviscid velocity distribution ^a^ nv ^ gc ^ ( j is calculated 
and throughout the whole report this is being referred to as V a *. The 
mass flow calculated at each station then is Q* and then Q* is compared 
with Q* real • 

The calculation of the density p* needed in the continuity equation 
is performed as follows: The absolute velocity is calculated as 



V * = M > u *) 2 + (V *) 2 + ( v r *) 2 



= /(V u *) 2 + (V a *) 2 + (V a *) 2 tan 2 ? 



( 19 ) 



The static temperature is calculated as 



The velocity of sound is calculated as 



( 20 ) 



* = _a_ = /YRgT 



( 21 ) 



The Mach number is calculated as 



M V* 

M = — 
a* 



( 22 ) 
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Then the density is: 



P* 



_P 

P atm 



C 

P 



R g H p atm 






1 

Tv^D 



n P * P / 

t atm / n Y-l 

■p p p l 1 ■*" p 

« H * w ^m Patm V 




1 

(Y-l) 




H* 

atm 




(23) 



The calculation of the mass flow rate at each radial position 
is performed with the subroutine SUMAT which uses Simpson's rule for 
unequal intervals. The subroutine SUMAT makes use of the subroutine 
INTP/) for interpolation which is described in detail in ref. 2. 

Once these calculations are performed the program gives control 
to the subroutine ITEEN which performs the iteration, until convergence 
has been realized and the continuity equation is satisfied. 

Block Q) 

The radial position of the streamlines at station (l) is considered 
the same throughout the calculation and in this step the radial positions 
of the streamlines at stations (2) and (3) are found with the condition 
that the same mass flow is allowed to pass through each streamtube with 
the newly calculated axial velocity distribution. 

Block Q 

The square root of the sum of the squares of the difference between 
the old and new position of the streamlines is calculated for stations (2) 
and (3). This is considered here as an indication of the convergence of 
the procedure. Additionally in this step all pertinent quantities that 
have not been calculated up to now are calculated (static and total pressures, 
relative velocities and mach numbers, angles, etc). 
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Block @ 



Once the prescribed number of iterations has been performed, 
some additional pertinent dimensional and non-dimensional quantities 
are calculated (a description of the calculation is given in Appendix C), 
and all the results are printed. 

The whole procedure is executed as many times as additional set 
of data exist. 

8. Description of the Use of the program 

This program has been constructed as a complement to ref. 2. An 
effort was made to take into account entropy and energy gradients (entropy 
gradients existing at the inlet or introduced through the loss correlations 
after each row, and aieigy gradients existing at the inlet or introduced by 
a non-uniform work distribution in the rotor) and compressibility effects. 

The non-dimensionsalization proposed in ref. 2 which does not take 
these effects into account is rendered thus incomplete and it was decided 
to carry out the calculation in the non-dimensional form proposed in 
reference 2 and modified slightly as described already, introducing, however, 
the data in dimensional form, reflecting thus the Mach number level. 

A detailed description of the program has already been given. The 
meaning of the weight factor has to be explained here. It happens sometimes 
in complicated cases that the iteration procedure diverges when the calculated 
corrections for the new streamline position are used in the whole. If, however, 
a fraction of the corrections is considered and introduced for the next itera- 
tion, the iteration procedure may be forced to converge. The weight factor 
( 0. < W.F. ^ 1.) introduced as data to the program defines the fraction of 
the correction to be used for the following iteration loop. 

Cases have been already run for no curvature effects and no losses. 

Then losses were introduced and curvature effects. It was found that to 
take into account the curvature effects and have an converging iteration 
process, a value of the weight factors smaller than unity ought to be used. 
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To facilitate the use of the program a table with a typical input 
has been prepared and given in Table II . The maximum number of stream- 
lines is taken to be eleven. The listing of the program along with the 
subroutines in use is given in Table III. The results of the already given 
typical input are given in Table IV. The explanation of the symbols 
used along with all the pertinent parameters used in the program are given 
in the FORTRAN Symbol Table. The output symbols not described in the 
FORTRAN Symbol Table will be found in the Output Symbol Table . For 
dimensional quantities the already given dimensions are used. For non- 
dimensional quantities the already described non-dimensionalization has 
been done. In the FORTRAN Symbol Table the dimensional variables are 
given. The non-dimensional ones are denoted in the program with the letter 
S at the end of the name of the dimensional quantity unless otherwise stated. 
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TABLE OF FORTRAN SYMBOLS USED IN THE MAIN PROGRAM 
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DHDRS(ll,3) t — -- total enthalpy gradient 
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atmospheric total enthalpy 



SYMBOL UNITS DESCRIPTION 
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R0(ll,3) P slug/ ft 3 static density table 
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Output Symbol table 
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APPENDIX A 



THE REACTION FACTOR WITH AXIAL VELOCITY VARIATIONS 




FIGURE 1A 



Rl £ R 2 
AH = oj(R 2 V u2 
APr ■ 



Rl v Ul ) 



r* = 



AP, 



ST 



Assuming no losses, for 
incompressible flow 



APr P 2 - Pi (Rtp^REL' ^tl^REL 



+ 



p ^ 2 2 - u : 







= i^Wu/.W^S + V al 2- Va2 2) 
+ i'-(U 2 2 -U 1 £ ) 



ap $t - »(r 2 v U2 . Rl v ui ) p 



which turns out to give a rather complicated formula for the reaction 
factor. Consequently, we shall assume that the quantity 



Up + U 2 

is specified along corresponding radii. 
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APPENDIX B 



The correlation described in reference 3 will be used for the 
calculation of losses (see also reference 4) 

The loss coefficient Q p is defined for rotors as 



C (APq)r 
p R - P tl - 



(APq*) R 

P , * _ P * 

^t! *1 



(Bl) 



where AP^ is the mass averaged total pressure loss, and for stators as 



C ( ap q) S ( af o*>S 

P S " p t 2 - p 2 ° p t 2 * - p 2 * 



This definition holds for incompressible flow where the density is assumed 
constant 

The diffusion factor D is given for rotors as 



W 2 R^ - Rg 

d r = 1 ' + a R (R x + R 2 ) w-l 

and for stators as 

v 3 R 2 ^ u 2 “ R 3 V u3 
D S = 1 - + dg (Rg + r 3 ) v 2 



(B3) 



(b4) 



or 



w* R*W, * _ P * W * 
n . W 2 *1 U 1 K 2 ^2 

d r - 1 " w * + 



U R (R-* + R 2 *j w x * 



(B3a) 



then 



p * v _ p -* v * 
z u 2 3 U3 

s * Vo* + (Ro* + R-3*) V * 



Vo* 

D =1 — 

V 2 ~ 



S v 2 3 7 2 



(B3b) 



r R 



jo.0O4 + 0.0639 (Dr + 0 . 1 ) 2 -» + 0.228D r 2 ' 02 

[1-Xr] 3 ' 77 ] (B5) 



where = 



R 2t " R 2 R2 t* ‘ R £ 



R R2 t ' R 2h R 2t* ‘ R 2h" 



(B 6 ) 
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For angles P 2 § 45°, see Fig. 1A, a correction is offered, where 
instead of Cpp the loss coefficient (^P R ) CQR is considered, where 



^PrW ” b P R cos p 2 /2 [l ~ H + IH§- 



(B7) 



For stators 



2a c 



= Ps 



co sa. 



|o.oo4 + 0.0639 (c s + o.i) 2 ^ 1 + 0.057 D s 2 ' 02 

[l- 1 s 1 3 - 7 ^J(b8) 



where no correction is offered for > ^5°. 

The shock losses are calculated as follows: 
The amount of supersonic turning is 

Au _ 9 ^ — - - (p _ p ) for the rotor 
R X 



(B9) 



0.625 



(og - <>3) for s i a i° r 



(BIO) 



The peak suction surface Mach number is then obtained from 

(M r ) = 1.095 + 0.03395 AU + 1.086 (M r - l.OO) 1 * 372 (Bll) 



for a rotor 



(Mg^) = 1.095 + 0.03395 Au + 1.086 (Mg - 1.00) 



1.372 



( B12) 



for a stator 

The shock losses are then calculated on the basis of the mean Mach number 



M 



M 



% + M s u 

2 

M 0 + Me 
2 b u 

2 



for a rotor 



for a stator 



(B13) 

(Bl4) 
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as 



S>H 



_ Ry 1 -!)^ 2 " |V('V- 1 ) r Y +1 " 1 

L(y-1)M 2 +2 J L2yM 2 -(y-1)J 

1 - [i^m 2 ]- y / ( ''- 1) 



1/(y-1) 



(B15) 



where the final loss coefficient is obtained as 



C = Op + c 



SH 



If the inlet Mach number is smaller than unity, the Mach number Mg u 
is calculated by taking the inlet Mach number M R ^ or M 2 to be unity. 

The mean Mach number M then is calculated using the actual inlet Mach 
number Mj^ or M^. If M < 1 no shock losses are assumed to exist. If M > 1. 
then the shock losses are calculated using equation (B15) and M = 1. 

As pointed out earlier, this procedure is adopted only if the inlet Mach 
numbers are smaller than unity. 
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APPENDIX C 



The following additional calculations are performed in the program, 
(a) The total pressure increase is 



-* * * 

AP = P - P 

t R t 2 t l 



* * * 

= P t " P t 
ST ^3 t l 



for the rotor 



for the stage 



(b) The static pressure increase is 

for the rotor 



* * * 
= P 2 - P 1 



* * * 
= p 3 ' p 2 



* * * 

AP = P - P 

ST 3 1 



for the stator 



for the stage 



(Cl) 

(C2) 

(C3) 

(c4) 

(C5) 



(c) The efficiency is calculated as follows 



Cp(T. 






Tl = 

's T 0 - T, 



- T ) 
t t ; 

2. t l 

is 

AH 
- T« 



total to total efficiency for the 
rotor (C6) 

static to static efficiency for the 
stator (C7) 



Cp(T, 



(^t) 



IS 



■v 



T-T 



(T 



(It) 



s-s 



AH 



- T ) 
3. 1 J 
is 

T - T 

3 1 



total to total efficiency for 
the stage (C8) 

static to static efficiency for 
the stage (C9) 



The calculation of the temperatures is being done as follows 

-(S -S )/Cp -(s*-s*) 

T t 0 = T t o e = T t o e 

2.2 2 

is 



(CIO) 
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(Cll) 



-(S -S )/Cp -(S*-S*) 

T = T e 0 = T e 0 ^ 

^is * ^ 



* *, 



r t = T t 6 

3- 3 

IS 



-(s 3 -s 1 )/c P -(s^) 

“ T, e 



(d) Mass averaged quantities are calculated. Assuming 

want to calculate the mass averaged value of the quantity Y ? 

R.. R* 

tip tip 

J 2tt P R Va Y dR J* 2Trp*R*V a * YdR* 

Y = ^hu b _ ^hub 



Q 



^real/S 



or 



. ^ 

^ * 

^real 



* 

R tip 

riub 



*T,* Tr * 



p R V a Y dR' 



* 



(C12) 

that we 
then 



(C13) 



(C13a) 
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SYMBOLS 



V absolute velocity 

W relative velocity 

U peripheral velocity 

H total enthalpy 

h static enthalpy 

Pj. total pressure 

P static pressure 

S entropy 

R radial distance from compressor axis 

w angular velocity 

L axial length (see Fig. 2) 

X angle between axial direction and streamline direction 

at stations (l), (2), and (3). 

6 R,AR radial distances defined in Fig 2. 

F,G functions defined in equations (l), (2), and (3) 

Cp specific heat at constant pressure 

z axial distance 

K constant used for the evaluation of curvature effects 

(see equation (4)). 

p static density 

p_k stagnation density 

a absolute angle measured from the axial direction 

(3 relative angle measured from the axial direction 

R gas constant 

g 

Y ratio of the specific heats at constant pressure and constant 
volume 
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T static temperature 

total temperature 
RF reaction factor 

AH total enthalpy increase in the rotor 

AP q loss in total pressure 

wall boundary layer blockage factor 
Q, mass flow rate 

a velocity of sound 

M Mach number 

AP static pressure increase 

AP^. total pressure increase 

D diffusion factor 

a solidity (chord to pitch ratio) 

Xr,X 3 non-dimensional quantities defined in equation B6 

T] efficiency 

Cp profile loss coefficient 

Cgk shock loss coefficient 

Au amount of supersonic turning 

M average mach number 

C total loss coefficient 

Subscripts 

1 station ahead of the rotor 

2 station behind the rotor and ahead of the stator 

3 station behind the stator 

a axial direction 

u peripheral direction 
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r radial direction 

is isentropic 

e entrance (see Fig. 

R rotor 

s stator 

h hub 

t tip 

m mean 

atm atmospheric 

ref reference 

st stage 

su suction surface 

cor corrected 

rel relative 



(start) 

Read required Data 

Set up a first approximation of the streamlines — (T) 

' i 

From distribution of R.F. and AH* calculate Vfl , Vu 2 ,Vu: 





Assume no losses in rotor and 
stator and no curvature 
effects 



Calculate Curvature terms, losses and entropy distribution ■ 


1 

1 


ir + 




Calculate entropy and enthalpy gradients and from them the 
functions F*and G 4 


-© 



3 

Solve radial equilibrium equation assuming a value for Vaj 



Test continuity —(7) 




Calculate new values of streamline positions (R* , R*3 ) —(8) 

Calculate error for positions (2) and (3 ) , print it and _/r\ 
calculate additional pertinent quantities. 




Calculate additional 
quantities and print — @ 
results 




STO F! 



NOTE* This procedure is repeated as many times as necessary to 
cover all cases introduced as data along with the program 



GENERAL FLOW DIAGRAM OF THE PROGRAM 
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